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ABSTRACT 


ww 


During  the  annual  Gray  Whale  migration  from  the 

Aleutians  to  Baja  California,  the  mammals  travel  in  coastal 

waters,  thereby  presenting  an  opportimity  for  the  study  of 

their  sound  spectral  and  source  levels  and  vocabulary.  How- 

< 

ever,  such  measurements  are  distorted  by  surface  and  bottom 
reverberation.  Using  the  theory  of  rough  surface  scattering, 
knowledge  of  the  bottom  impedance,  and  correlation  techniques, 
it  is  possible  to  decompose  the  shallow  water  reverberation 
into  the  contributions  from  different  paths . From  this , the 
range,  depth  and  the  deverberated  spectral  source  levels  of 
the  sounds  of  the  mammal  can  be  determined  by  use  of  only  one 
hydrophone  rather  than  the  conventional  three  or  four.  The 
theory,  deverberation  programming,  and  experimental  results 
are  presented  for  a model  of  the  whale's  pulsed  radiation  in 
a laboratory  model  coastal  environment. 
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I . INTRODUCTION 


Once  yearly,  the  Gray  Whale,  Eschrichtius  glaucus , 
migrates  southward  from  the  Aleutians  and  passes  very  close 
to  the  California  coast  in  shallow  water.  During  this 
migration,  the  sounds  in  the  water  close  to  the  whales  can 
easily  be  recorded;  but,  they  may  not  be  the  true  sounds 
produced  by  the  whale.  The  whale  normally  produces  an 
intermittent,  short  duration  signal  which  in  shallow  water 
is  received  at  the  hydrophone  via  direct,  surface  reflected 
and  bottom  reflected  paths.  Since  the  path  lengths  are 
different,  the  signals  arrive  at  the  hydrophone  at  different 
times  and  they  interfere.  To  obtain  the  true  somds  produced 
by  the  whales,  this  interference,  which  is  called  reverbera- 
tion, must  be  eliminated. 

The  purpose  of  this  thesis  is  to  model  in  the  laboratory 
the  whale  sounds  in  the  shallow  water  and  to  develop  a 
technique  to  determine  the  range,  eliminate  the  surface  and 
bottom  reverberation  and  calculate  the  spectral  source  levels 
as  a function  of  time. 


II.  THEORY 


In  a reverberant  environment,  the  original  signal  can 
only  be  realized  if  the  reverberation  is  removed.  The 
method  used  to  remove  the  reverberation,  which  is  here 
called  "deverberation, " assumes  that  the  whale  is  a point 
source,  the  geometrical  spreading  is  spherical,  the  water  is 
isovelocity  and  the  water  depth  is  constant. 

Before  the  deverberation  technique  can  be  applied,  the 
direct,  surface  reflected  and  bottom  reflected  path  dis- 
tances must  be  known.  Normally,  this  information  is  obtained 
by  knowing  the  source  position  and  calculating  the  path 
distances  from  the  known  geometry.  Determining  the  hori- 
zontal source  position  is  difficult  and  requires  at  least 
two  directional  or  three  omnidirectional  hydrophones,  accu- 
rately fixed  with  respect  to  each  other  at  all  times.  To 
determine  the  depth  requires  at  least  one  additional  hydro- 
phone at  a different  depth.  Generally,  three  or  four 
hydrophones  are  deployed  CRefs . 1 and  2].  With  each  addi- 
tional hydrophone,  the  complexity  of  the  system  is  increased 
since  each  hydrophone  requires  its  own  amplifying  and  record- 
ing system.  In  shallow  water,  however,  taking  advantage  of 
the  surface  and  bottom  reflections,  one  can  use  a single 
hydrophone  and  gather  all  the  information  necessary  for  the 
application  of  the  deverberation  technique. 

Consider  the  direct,  surface  scattered,  and  bottom 
scattered  sounds  received  at  only  one  hydrophone  which 
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is  deployed  in  shallow  water  as  depicted  in  figure  1.  It 
will  be  shown  that  when  the  differences  between  the  arrival 
times  for  the  three  paths  are  known  the  three  path  distances 
can  be  calculated.  Using  the  surface  specularly  scattered 
path,  R it  is  seen  that 

O 

t 11 

= D’  + s'  ^ 

where 

s’*  = 7— [r*  - (.H-dV 

s'‘^  = [r*-  - (.H-D)' 

CH  + D)‘  ^ 

Substituting  and  rearranging  terms  gives 

(r*.  4HD)^ 

Using  Tg,  the  time  difference  between  the  direct  path 
arrival  and  the  surface  reflected  path  arrival  and  C,  the 
mean  speed  of  sound,  gives 

C Ts  = R’s  + R s - R 

and  therefore 

X 

(3)  CT^  + 4V\d)^  - R 

Similarly  using  the  bottom  specularly  scattered  path 

rb  - r'b  ^ r; 

R'a  = Ci-D')*  + B'* 

-D)^  + B"* 

o 
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LOYMENT  OF  ONE  HYDROPHONE  IN 
SHMLOW  W^TER  ENVIRONMENT 


where 

I* 

■(&'»%>  I"'  ■ ■“■’’I* 

Substituting  and  rearranging 

T!'b  °'(A^-W-*Py  + 

I 

.5.  o'*  [^^4  4 (z^+ WD  - -id)]  ~ 

(5)  ^21  -H  -D)  i 

Using  Tg,  the  time  difference  between  the  direct  path 
arrival  and  the  bottom  reflected  path  arrival  gives 

and,  therefore, 

X 

(6)  cT«  = \'R^  + + - R 


Solving  equations  (3)  and  (6)  simultaneously  for  R,  the 
range  of  the  source  from  the  receiver  and  D,  the  depth  of 
the  source,  produces 

,,,  r^-  -v  4ZH)  - T^CcTTb)^. 

<i-Lhtb + rsU-H)] 

fa^  D " 4HD  ~ 

(8)  R - 2CTs 

The  equation  for  the  range  is  left  as  a function  of  the  water 
depth  to  facilitate  its  being  programmed.  Now  that  the  range 
and  depth  are  known,  the  other  two  path  distances  can  be 
calculated.  From  equations  (1)  and  (2) , the  surface  reflected 
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path  distance  is 


(9)  -V 

and  from  equations  (4)  and  (5) , the  bottom  reflected  path 
distance  is 

(10)  UD  - ZH  - ^D')J  ^ 

Therefore,  when  the  time  arrival  differences  for  the  differ- 
ent paths  are  known,  equations  (8) , (9) , and  (10)  can  be 
used  to  determine  the  first  three  path  distances.  The  paths 
for  multiple  reflections  can  be  calculated  directly  from  the 
known  geometry,  assuming  specular  scatter. 

For  a transient  signal,  determination  of  the  differential 
arrival  times  Xg  and  Xg  can  be  realized  by  the  use  of  an 
autocorrelation  technique.  The  correlation  function  can  be 
defined  as  [Ref,  3] 

<ii)  R (.t)  --  -ulI} 

where  v is  the  time  average,  u is  the  mean,  and  E is  the 
expected  value  of  the  received  signal.  This  process  is 
programmed  using  a digital  summation 

(12)  R Ct)  - j Wt  W ■ '^11. 

with  n representing  the  total  number  of  samples  in  the  record. 
Performing  the  autocorrelation  on  the  reverberant  signal  at 
the  one  hydrophone  gives  peaks  at  delay  times  corresponding 
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to  zero  delay  time,  and  the  arrival  delay  times  from  the 
reflected  signals.  The  peaks  are  realized  only  when  the 
direct  signal  is  delayed  enough  to  correlate  with  the 
reflected  signals.  The  computer  program  called  AUTOPEAK 
performs  the  autocorrelation  and  then  applies  an  envelope 
detection  to  determine  the  delay  times  for  the  peaks.  These 
delay  times  are  then  used  in  equations  (7),  (8),  (9),  and 
(10)  to  determine  the  range,  depth,  and  reflected  path  dis- 
tances. Thereby,  all  the  geometrical  information  necessary 
for  the  deverberation  technique  has  been  obtained. 

Computer  programs  have  been  developed  for  deverberation 
in  either  the  frequency  or  time  domain. 

The  computer  program  designed  to  eliminate  the  reverbera- 
tions in  the  frequency  domain  is  called  DEVERB.  For  the  time 
being,  assume  only  one  frequency  component,  whose  amplitude 
and  frequency  are  fxinctions  of  time.  For  a signal  with  time- 
varying  frequency  and  amplitude,  we  can  describe  the  pressure 
at  the  receiver,  due  only  to  the  direct  path  signal  as  [Ref. 
4] 

j uj  W i 


(13)  Pjj  (O  - c Ci')  e 


Then,  taking  into  account  spherical  divergence,  the  spatial 
phase  shift,  and  specular  scattering  from  a Gaussian  rough 
surface,  the  pressure  at  the  hydrophone  due  to  the  surface 
scattered  signal  can  be  written  as 


( 14)  = Pj  (j.-T.')  = f e * t W 

t 
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and  the  pressure  due  to  the  once  bottom  reflected  signal  is 

(15)  Pb  f ® t w e. 

■1^'Cb  ' 

and  e ^ are  the  frequency-dependent  pressure  amplitude 
reflection  coefficients,  and  V and  tt  are  the  phase  shifts 
due  to  the  bottom  and  surface  reflections,  respectively. 

The  surface  reflection  coefficient  depends  on  the  roughness 
of  the  surface.  The  exponent  for  specular  scattering  is 
[Ref.  5] 


(16)  9 


z - 


4Tr  cr 


~ cos  ©s 


where  CT  is  the  R.M.S.  wave  height,  X is  the  wavelength  of 
the  sound  signal  and  0g  is  the  angle  of  incidence.  Equations 
(14)  and  (15)  are  derived  from  previously  received  direct 
path  pressures,  corrected  for  path  differences  and  phase 
shifts  and  represent  the  pressures  due  to  the  single  reflected 
paths.  The  coherent  sum  of  equations  (13),  (14)  and  (15)  is 
the  pressure  sensed  by  the  receiver  when  these  three  compo- 
nents are  present. 

The  signal  processing  is  done  digitally.  Therefore,  the 
continuous  time  dependence  is  replaced,  whenever  it  occurs 
in  the  previous  equations,  by  digital  block  numbers,  indi- 
cated by  the  subscript  index  "k.”  Each  block  contains 
enough  data  samples  of  the  incoming  signal  to  give  the  desired 
spectral  frequency  resolution  during  a block  duration  vhich 
is  small  compared  to  the  total  duration  of  the  time-varying 
signal.  The  frequency  change  of  the  signal  within  a block 
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duration  is  assumed  to  be  much  smaller  than  the  frequency 
resolution  of  the  digital  processing.  The  block  duration  is 


T,  ^ NUMBER  OF  SAIiPLES  IN  THE  BLOCK 
SAMPLING  FREQUENCY 

The  time  the  block  ends  is  related  to  the  continuous  time  t 


(17)  t = TK;  K = 1,  2,  3 • • • 

The  index  "i"  is  used  for  the  spectral  frequency  compo- 
nent of  the  complex  wave  being  analyzed. 

The  equation  for  the  source  pressure  at  unit  distance 
using  the  frequency  deverberation  correction  is 

- tTC.} 


where  l(K-N),  l(K-M)  and  l(K-L)  are  unity  factors  with  values 
l(K-N)  =1  K > N 

= 0 otherwise 

l(K-M)  =1  K > M 

= 0 otherwise 


l(K-L)  = 1 
= 0 


K > L 


otherwise 


ttl 

The  pressure  amplitude  of  the  i*"  frequency  component  in 
block  K of  the  receiver  reverberant  signal  is  represented 
by  ^ and  its  phase  represented  by  4)^^  The  deverberated 
pressure  amplitude  is  represented  by  Dj^  ^ and  its  phase  by 
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The  first  term  on  the  right  hand  side  of  the  equation 


represents  the  received  signal,  The  second  term  represents 
the  correction  due  to  a single  specular  scatter  from  the 
surface;  the  third  represents  a single  bottom  reflection 
correction,  and  the  fourth  represents  the  correction  for  a 
path  which  includes  one  surface  and  one  bottom  reflection. 

The  equation  can  be  expanded  to  include  other  multiple  reflec- 
tions . 

The  block  indices  are  determined  by 


(19)  M = K - ^ 

T 

(20)  N = K - 

T 

(21)  L = K - 

T 

The  output  of  the  frequency  deverberation  program  is  a 
series  of  spectra  of  the  consecutive  blocks  and  its  Fourier 
transform  is  a time  plot  of  the  deverberated  signal. 

The  above  procedure  takes  the  signal  from  the  time 
domain  (the  time  series  after  A/D  conversion)  by  Fourier 
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transform  to  the  frequency  domain,  where  the  known  frequency 
dependent  reflection  coefficients  are  easily  applied,  and 
then  back  to  the  time  domain  to  verify  the  effectiveness  of 
the  process. 

When  the  reflection  coefficients  can  be  assumed  to  be 
frequency-independent,  a simple  point-by-point  deverberation 
procedure  can  be  applied  in  the  time  domain.  The  applicable 
temporal  deverberation  equation  is 

represents  the  pressure  amplitude  for  the  sample.  The 

other  terms  are  similar  to  those  in  equation  (18) . For  low 
roughness  surfaces  g < 1 and  the  use  of  over  the 

appropriate  frequencies  will  be  a good  approximation  which 
is  essentially  independent  of  frequency.  One  advantage  of 
the  temporal  deverberation  technique  is  the  relative  freedom 
from  restrictions  of  block  size;  the  block  size  is  determined 
only  by  the  desired  frequency  resolution  and  the  rate  of 
change  of  frequency  of  the  transient  source. 
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III.  PROCEDURE 


In  order  to  model  the  Gray  Whale's  environment  in 
Monterey  Bay,  a three  meter  cube  "anechoic"  water-filled 
tank  was  used.  An  artifical  bottom  made  of  hard  rubber  with 
an  experimentally  determined  pc  of  approximately  2.4  x 10 
mks  rayls  was  positioned  one  meter  below  the  water  surface. 
This  type  of  bottom  was  chosen  for  its  specific  acoustic 
impedance  since  the  bottom  at  the  listening  area  in  Monterey 
Bay  has  a pc  of  approximately  3 x lO'^  mks  rayls.  The  depth 
of  the  bottom  and  relative  placement  of  the  source  and 
receiver  were  determined  in  order  to  obtain  realistic  modeled 
delay  times  between  the  reflected  signals  and  the  direct 
signal.  A 1.8  cm  diameter  spherical  source  was  used  because 
of  its  small  size  and  its  ability  to  transmit  a signal  with 
minimum  distortion;  but  it  also  limited  the  minimum  frequency 
to  about  10  kHz.  An  FM  up-sweep  of  varying  widths  and  sweep 
rates  was  used  to  model  one  of  the  sounds  produced  by  the 
Gray  Whale. 

The  equipment  was  connected  as  in  figure  2 with  the  mas- 
ter unit  pulse  generator  being  used  to  synchronize  the  samp- 
ling frequency  oscillator  which  determined  the  start  and  stop 
frequencies,  sweep  rate  and  pulse  width.  The  pulse  repeti- 
tion rate  was  determined  by  the  master  unit  pulse  generator. 
The  signal  was  then  amplified  and  sent  to  the  source.  The 
reverberent  signal  was  received  by  an  LC-10  hydrophone,  and 
then  amplified  again  and  sent  to  the  analog- to-digital 


FIGURE  2:  EQU'PtAENT  BLOCK  DIAGRKM 


converter.  The  A/D  converter  was  triggered  through  the 
delay  generator  via  the  unit  pulse  generator  and  the  sampling 
frequency  oscillator.  A sampling  frequency  of  320  kHz  was 
used,  The  delay  generator  delays  the  trigger  by  the  time 
required  for  the  signal  to  be  transmitted  through  the  water 
and  then  this  delayed  pulse  triggers  the  unit  pulse  generator. 
The  unit  pulse  generator  was  used  to  determine  the  total  on- 
time  of  the  sampling  frequency  oscillator.  The  oscillator 
determined  the  samples  per  second  taken  by  the  converter 
during  the  oscillator's  on-time.  This  complex  equipment  set- 
up was  designed  to  allow  the  A/D  converter  to  sample  only 
during  the  time  when  the  direct  and  reflected  signals  were 
present,  thereby  reducing  the  amount  of  computer  time  and 
storage  required  to  process  the  data.  Reflections  from 
other  surfaces  in  the  tank  were  eliminated  wherever  possible 
in  the  sampling  time  by  varying  the  pulse  repetition  rate 
of  the  master  pulse  generator  or  by  varying  the  source  and 
receiver  placement.  After  the  analog  signal  was  changed  by 
the  converter  to  digital  form,  it  was  stored  on  cassette 
memory  and  then  analyzed  using  the  programs  AUTOPEAK  and 
DEVERB . 


IV.  DATA  PROCESSING  AND  RESULTS 
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The  autocorrelation  plots  of  reverberation  for  two 
different  sweep  widths  are  seen  in  figure  3 with  the  scale 
factor  for  the  delay  time  equal  to  3.125  y seconds.  The 
plot  of  the  90  kHz  sweep  width  shows  a steeper  slope  of  the 
envelope,  thus  a more  clearly  defined  peak  than  the  one  with 
only  a 10  kHz  sweep  width.  This  indicates,  as  expected,  that 
as  the  difference  between  the  upper  and  lower  frequencies 
decreases,  the  correlation  peaks  become  harder  to  determine. 

In  the  limit  of  only  one  frequency  being  present,  there  would 
be  no  peaks  in  the  envelope.  Figure  4 shows  the  range  error 
(computer  determination  compared  to  direct  measurement)  versus 
the  ratio  of  the  upper  frequency  to  the  lower  one.  The  graph 
indicates  that  for  a ratio  above  1.2  to  1 the  frequency  sweep 
of  the  signal  is  sufficient  to  get  accurate  time  differences 
for  the  reflections  and  thereby  to  determine  the  source  range 
and  depth  from  the  autocorrelation  processes. 

The  frequency  deverberation  program  is  designed  to  give 
a true  spectrum  of  the  signal  by  eliminating  frequency  depen- 
dent reverberations.  Since  the  signal  is  time-varying,  a 
small  block  of  time  is  desirable  to  keep  the  change  in  the 
signal  to  a minimum  diaring  the  block.  The  limiting  factor 
to  the  minimum  block  size  is  the  desired  frequency  resolu- 
tion, For  a spectrum  to  accurately  represent  an  instant  of 
time  rather  than  being  a spectrum  averaged  over  a length  of 
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FIGURE  3. 

AUTOCORRELATION  OF  DIRECT  AND  SURFACE  REFLECTED  CHIRP 
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time,  the  frequency  change  during  the  block  time  should  be 
much  less  than  the  frequency  resolution.  In  addition,  where 
possible,  the  block  duration  should  be  equal  to,  or  a sub- 
multiple of  I - T^l  in  order  to  permit  equation  (18)  to 
be  applied. 

From  figure  4,  it  is  known  that  depth  and  range  can  be 
determined  only  when  the  frequency  ratio  of  the  chirp  is 
greater  than  1.2  to  1,  Figure  5 shows  the  reverberant  and 
deverberated  signals  for  a chirp  from  46  kHz  to  54  kHz  with 
a ratio  of  1,17  to  1,  The  reverberant  signal,  top  of  figure 
5,  was  divided  into  blocks,  equal  to  Tg  - Tg,  and  then  trans- 
formed back  into  the  time  domain  which  is  shown  at  the  bottom 
of  the  figure.  Some  improvement  can  be  seen  but  the  expected 
slowly  increasing  frequency  at  the  approximately  constant 
amplitude  has  not  been  realized.  It  is  believed  that  the 
inadequate  deverberation  is  due  to  the  fact  that  the  frequency 
sweep  during  a block  is  approximately  four- tenths  of  the 
frequency  resolution.  A slower  sweep  rate  or  a larger  block 
duration  would  have  cured  this  problem.  However,  a slower 
sweep  rate  for  the  model  would  have  decreased  the  accuracy  of 
the  range  determined  by  the  autocorrelation;  and  a Targe  block 
duration  was  precluded  by  the  geometry  which  determined 
I Tg  - ^^3!  . The  temporal  deverberation  technique  which  is 
presented  next  did  not  suffer  from  these  limitations. 

The  result  of  applying  the  temporal  deverberation  program 
to  a signal  with  a sweep  width  from  5 kHz  to  95  kHz  can  be 
seen  in  figure  6 with  the  reverberant  signal  on  the  top  and 
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the  deverberant  signal  below,  Since  the  source  was  resonant 
at  64  kHz,  the  FM  sweep  grows  in  amplitude  to  64  kHz  and 
then  decreases.  The  deverberated  signal  shows  the  frequency 
and  amplitude  modulation  cleanly  until  the  end  of  the  direct 
signal.  The  reverberation  after  that  time  is  due  to  scatter- 
ing from  the  side  walls  of  the  tank.  Possible  reverberation 
due  to  the  fourth  and  later  terms  on  the  right  hand  side 
of  equation  (22)  were  excluded  by  limiting  the  duration  of 
sampling  by  the  computer. 
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V.  CONCLUSIONS 


The  work  described  in  this  thesis  demonstrates  the 
feasibility  of  obtaining  a non-reverberant  spectrum  of  a 
transient  source  in  a reverberant  environment.  The  technique 
includes  calculation  of  the  autocorrelation  of  the  received 
signal  to  determine  range  and  depth  of  the  source  and  com- 
puter processing  to  correct  for  the  surface  and  bottom  reflec- 
tions . 

The  autocorrelation  function  provides  an  accurate  method 
for  obtaining  the  range  and  depth  of  a source  of  transient 
sound  in  shallow  water.  The  correlation  technique  can  be 
performed  for  a chirp  sound  with  ratio  of  upper  to  lower 
frequency  of  greater  than  1.2  to  1.  At  least  two  reflections 
are  required  to  obtain  the  depth  and  range  of  the  source  with 
respect  to  the  receiver. 

Frequency  and  time  deverberation  programs  which  use  the 
position  data  from  the  autocorrelation  have  been  developed 
to  eliminate  the  reverberations  and,  thereby,  to  obtain 
corrected  spectra  or  corrected  time  plots . Either  technique 
can  be  used;  however,  because  the  output  of  the  computer  is 
a time  series,  it  is  natural  to  apply  temporal  deverberation. 
This  becomes  very  simple  if  the  surface  and  bottom  reflection 
coefficients  are  independent  of  frequency. 
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REM 
REf- 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
DIM 
DIM 
DIM 
D I M 
REM 
Z4$="N 
Z$=«Y" 
REM 
REM 
GOSUB 
IF  Z4$ 
REM 
GOSUB 
REM 
GOSUB 
REM 
GOSUB 
REM 
GOSUB 
IF  Z3$ 
REM 
GOSUB 
REM 
GOSUB 
IF  Z8= 
REM 
GOSUB 
REM 
GOSUB 
REM 
GOSUB 
REM 
GOSUB 
REM 
GOSUB 
PRINT 
INPUT 
IF  Z$= 
PRINT 
INPUT 
PRINT 
INPUT 
IF  Z$= 
Z2$="Y 
GOTO  1 
REM 
REM 
IF  Z$ 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
INPUT 
PRINT 
PRINT 
PRINT 
INPUT 
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AITOPEAK  — 9/29/77 

JEANIE  SAVAGE,  PROGRAMMER 
SPECIFICATIONS  BY  RICK  ECSTIAN 
TFIS  PROGRAM  PERFORMS  AN  AUTO- 
CORRELATION ON  SIGNAL  DATA  AND 
THEN  PICKS  OUT  THE  TWO  PEAK 
VALUES. 


(48)  ,V1$(6)  , Z$( 1) 

1000) ,6(1000) ,Y( 1000) 

S(  1),Z2$(1),S(1),D(1),Z4$(1) 
500,1) ,Z3$(1),M (50,2),  R(  1,2) 


=!>**DRIVER  ROUTINE*’*'* 

PERFORM  INITIALIZATION  PROCEDURE 

500 

=«Y"  GCTO  140 

READ  DATA  FROM  TAPE 


1000 

1300 

1500 


BUILD  CTHEP  ARRAY 
PERFORM  AUTO-CORRELAT ICN 
DETERMINE  INTERMEDIATE  PEAK  VALUES 


2000 
= ''N>'  GOTO  185 

PRINT  INTERMEDIATE  VALUES 

2500 

DETERMINE  TWO  PEAK  VALLES 

3000 

1 GOTO  150 

PRINT  TIME  DIFFERENCES 


4000 

4500 

5000 

5500 


CALCULATE  SOURCE  DEPTH,  S-R  RANGE 
PRINT  DEPTH  AND  RANGE 
CALCULATE  REFLECTION  PATH  DISTANCES 
PRINT  DISTANCES 


6000 

"ARE  YCU  FINISHED?  (Y  OR  N) — " 

Z$ 

"Y"  THEN  STOP 

"SAME  DATA?  (Y  OR  N) — " 

Z4$ 

"SAME  PARAMETERS?  (Y  OR  N) — " 

Z$ 

"N"  GOTO  115 

II 
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***INI TIALIZATION  PRCCECURE^^^ 

"N"  GOTO  525 
"AUTOPEAK" 

"SAMPLING  FREQUENCY  MUST  BE  FOUR  TIMES  THE" 

" GREATEST  FREQUENCY  OF  INTEREST." 

"NUMBER  OF  POINTS  PER  BLOCK--" 

N2 

"MINIMUM  TIME  DIFFERENCE  BETWEEN  DIRECT  PATH" 
" AND  FIRST  REFLECTED  PATH.  RECEPTION  (IN" 
PRINT  " SAMPLES) — " 

D8 
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535 
540 
545 
55C 
555 
560 
56  2 
562 
565 
57C 
575 
58C 
585 
590 
595 
60C 
605 
61C 
61  5 
62C 
625 
630 
1000 

1004 

1005 
1010 
1015 
1020 
1025 
1030 
1032 
1035 
1040 
1045 
105C 
1055 
1060 
1065 
107C 

loac 

1085 

1090 

1095 

1100 

1300 

1305 

1310 

1315 

1320 

1325 

1500 

1505 

1510 

1515 

1520 

1525 

1530 

1535 

1540 

1545 

1550 

1555 

1560 

1565 

1570 

1575 

1580 

1585 

1590 

1595 

2000 

2005 


PRINT  "NUMBER  CF  POINTS  TO  BE  USED  FROM  EACH  PIOCK" 
PRINT  " (NUMBER  CF  POINTS  PLUS  DELAY  MUST  EE  LESS" 
PRINT  " THAN  NUMBER  CF  POINTS  IN  BLOCK) — " 

INPUT  N1 

PRINT  "MAXIMUM  LAG — " 

INPUT  LI 

PRINT  "PRINT  PEAK  VALUES?  (Y  OR  N)--" 

INPUT  23$ 

IF  Z$="N"  THEN  PRINT  "CONTINUING  V^ITH  CALCULAT ICNS" 

IF  Z$="N"  THEN  RETURN 

PRINT  "IS  THIS  THE  1ST  BLOCK  OF  A MULTIPLE  RUN?  (Y/N) — 
INPUT  22$ 

PRINT  "SAMPLING  FREQUENCY — " 

INPUT  SI 

PRINT  "SPEED  OF  SOUND  (M/SEC) — " 

INPUT  C 

PRINT  "DEPTH  OF  WATER.  (Mi  — " 

INPUT  H 

PRINT  "DEPTH  OF  RECEIVER  (M)--" 

INPUT  R 

RETURN 

REM 

REM  **’!'ROUTINE  TO  READ  DATA  FROM  TAPE*^^* 

PRINT  "SWITCH  TO  HIGH  SPEED." 

IF  Z2$="Y"  THEN  INPUT  ON  (2)  V$ 

K1=0 

FOR  1 = 1 TO  (N2/8) 

INPUT  ON  (2)  V$ 

FOR  J=1  TO  48  STEP  6 
Vl$="  •• 

Vl${ 1,6) =»" 

FOR  K=0  TO  5 

IF  V$( J+K, J+K)="  " GOTO  1050 
V1$=V1$+V$( J+K, J+K) 

NEXT  K 

A (K1)  = VAL{ Vl$) 

K1=K1+1 
NEXT  J 
NEXT  I 

PRINT  "CONTINUE?" 

INPUT  Zl$ 

PRINT  "CONTINUING  WITH  CALCULATI CNS 

RETURN 

REN 

REM  *«#BUILD  CTHER  ARRAY^^-*^ 

FOR  1=0  TO  Nl-1 
B(I)=A( I+D8) 

NEXT  I 
RETURN 
REM 

REM  *«>;«CROSS-CCRRELATiaN  ROUTINE*^* 

A8=0 
B 8 = C 

FOR  1=0  TO  Nl-1 
A8  = A8+A(  1) 

B8=88+B( I) 

NEXT  I 
A8=A8/N1 
B8  = B8/M 
FOR  1=0  TO  LI 
N9=N1-I 

se=c 

FOR  J=0  TO  N9-1 
I9=J+I 

S8=S8+(A(J)-A8)*(B(I9)-B8) 

NEXT  J 
Y( 1)=S8/N9 
NEXT  I 
RETURN 
REM 

REM  ***OETERMINE  INTERMEDIATE  PEAK  VALUES*** 

IF  Y(1)>Y(0)  THEN  E=l 
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2010  IF  Y(1J<Y(0)  THEN  E=-l 

2015  P9=-l 

202C  FOR  1=2  TO  LI 

2025  IF  e=(-l  ) GOTO  2060 

2030  IF  Y(1J>Y(I-1)  GOTO  2085 

2035  PS=P‘5+1 

20AC  P (P9,0)= I-1+C8 

2045  P tP9,l ) = Y( I-l) 

205C  E=-E 
2055  GOTO  2085 

2060  IF  Y(II<Y<I-1}  GOTO  2085 

2065  P9=P9+1 

2070  P (P9,0)=I-1+C8 

2075  P (P9,1)=Y( I-l) 

2080  E=-E 
2085  NEXT  I 
2090  RETURN 
2095  RE)< 

2500  REM  ***PRINT  INTERMEDIATE  PEAK.  VALUES*’!-’* 

2505  PRINT 

2510  PRINT  "“PEAK  VALUES — " 

2515  PRINT 

2520  FOR  1=0  TO  P9 

2525  IF  I/5=INT(I/5)  GOTO  2540 

2530  PRINT  P( 1, 1)  , 

2535  GOTO  2550 
2540  PRINT  P(  1,1) 

2545  PRINT 
2550  NEXT  I 
2555  PRINT 
2560  RETURN 
2565  REM 

3000  REM  ’!'*’KOET ERMINE  TWO  PEAK  VALUES’!’** 

3003  Z8=0 
3005  E=C 
3010  Kl=-1 

3015  FOR  1=0  TO  P9-3 
3020  IF  E=1  GOTO  3040 

3025  IF  A6S(P(I,1))>ABS(P(I+1,1))  OR  ABS(P(I,1))> 

ABS( P( 1+2,1 ) ) GOTO  3075 

3030  IF  A6S<P(I,1))>  ABS(P( 1+3,1) ) GOTC  3075 
3035  E=1 

3040  IF  ABS(P  ( 1 ,1  ) XABS  ( P(  1+1,1)  OR  AeS(P(I,l))< 
ABS(P( 1+2,1) ) GOTC  3075 
3045  IF  ABSlPt 1,1  ))<ABS(P( 1+3,1))  GOTO  2C75 
305C  E=C 
3055  K1=K1+1 
3060  M(K1,0)  = P{  1 ,0)/Sl 
3065  M{K1,1)=P(I,1) 

3070  M(K1,2)=I 

3075  NEXT  K 

3076  IF  K1>=0  GOTC  3081 

3077  PRINT  "CHANGE  LAG  — MAXIMUM  LAG=" 

3078  INPUT  LI 

3079  ze=i 
308C  RETURN 

3031  REM  DETERMINE  TWO  HIGHEST  PEAKS 
30E2  M8=M{C,1) 

3085  I9=C 
3090  FOR  1 = 1, K1 

3095  IF  ABS<M(I,1))>ABS(M8)  THEN  19=1 
3100  IF  A8S(M<I,1))>  AaS(M8J  THEN  MS=M(I,1) 

3105  NEXT  I 
3110  R(0,0)=M8 
3115  R(C,1)=M{ 19,0) 

3120  R(C,2)=M(I9, 2) 

3125  ie=c 

3130  IF  19=0  THEN  18=1 
3135  M6=M(I8, 1) 

3137  IF  Kl=i  GOTO  3165 
3140  FOR  1=18+1  TC  K1 
3145  IF  1=19  GOTO  3160 
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315C  IF  ABS(M(I,1))>A6S(M8)  THEN  18=1 
3155  IF  ABS(M(Itl))>ABS(M8)  THEN  MS=M(I,1) 

316C  NEXT  I 
3165  R(1,0)=N8 
317C  R(1,1)  = M(  18,0 
3175  R(1,2)=M( 18,2) 

3iaC  RETURN 
3185  REM 

4000  REM  ♦♦♦PRINT  TIME  01 FFERENCES**^ 

4005  PRINT 
4010  PRINT 

4015  PRINT  "TIME  DIFFERENCE  BETWEEN  DIRECT  AND  SURFACE 
REFLECTED  PATHS” 

4020  PRINT  " 

4025  PRINT  TAB(9)  ;''PEAK  VALUE” 

4027  A3  = R(0,1  )^10C0 

4030  PRINT  USING  •'  22,25)0)33  MSEC'SAi 

4035  PRINT 
404C  PRINT 

4045  PRINT  "TIME  DIFFERENCE  BETWEEN  DIRECT  AND  BCTTCM 
REFLECTED  PATHS" 

4050  PRINT  " 


4055  PRINT  TAE(9);"PEAK  VALUE” 

4057  A3=R(1 ,1 )*10C0 

406C  PRINT  USING  " 23.23333  MSEC”,A2 

4065  PRINT 

407C  PRINT 

4075  RETURN 

4080  REM 

450C  REM  ♦♦♦CALCULATE  SOURCE  DEPTH  AND  S-F  RANGE^-^^ 

4505  FOR  1=1  TO  1 

451C  GC=R(0,I)^(C^R(1,I))  2 

4515  G1=C  2*RiO, I )*R{ 1 , I)+4^H  2-4^H^R 

452C  G2=RI 0, I )«G1 

45  25  G3=4^(R(  0, I ) ♦ ( R-H ) -R*R ( 1 , 1 ) j 

453C  S( I-l)=< G0-G2)/G5 

453  5 C( I-1)  = 2^R^S(I-! ) / ( C«R ( 0 , 1 ) ) -C^R ( C , I ) /2 
454C  NEXT  I 
4545  RETURN 
4550  REM 

5C0C  REM,  ♦♦♦PRINT  DEPTH  AND  DISTANCE  VALUES^^^ 

5005  PRINT  "RANGE  OF  SOURCE  FROM  RECEIVER  AND  SOURCE  DEPTH 
IN  METERS:" 

5010  PRINT  " 


5015  PRINT  TAB(9};"PEAK  VALUES" 

5020  PRINT  TAB(5) ;”DEPTH”;TAB( 19) ;"RANGE" 

5030  PRINT  USING  " 33.23323  ”,S(0),D(0) 

5035  PRINT 
504C  PRINT 
5045  PRINT 
505C  RETURN 
5055  REM 

5500  REM  ♦♦♦CALCULATE  DISTANCES  CF  REFLECTED  PATHS** 

5505  REM  CALCULATE  TRANSVERSE  SOURCE-RECEIVER  DISTANCE 
5510  T=SQR(D(0)  2-(R-S(C))  2) 

5515  REM  CALCULATE  DISTANCE  OF  SURFACE  REFLECTION  PATH 
552C  X=S(0)*T/(R  + S(0)  ) 

5525  Y=F*T/ (R+S(0)) 

5530  U=SQR(X  2+S(0)  2) 

5535  W=SQR(Y  2+R  2) 

5540  DC=U+W 

5545  REM  CALCULATE  DISTANCE  OF  BOTTOM  REFLECTICN  PATH 
5550  A=(H-S(0))*T/(2*H-R-S(C)) 

5555  B=(H-R)*T/(2*H-R-S(0)  ) 

5560  E = SQR(A  2-*-(H-S(0))  2) 

55o5  F=SQR(B  2+(H-R)  2) 

5570  D1=E+F 
5575  RETURN 
558C  REM 


32 


-I 


>1 


6000 

RE 

M 

6005 

PR 

INT 

6010 

PR 

INT 

6015 

PR 

INT 

6020 

PR 
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603C 

PR 
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6035 

PR 
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6040 

PR 

INT 

6045 

PR 

INT 

6050 

PR 

INT 

6055 

RE 

TURN 

6060 

RE 

M 

6065 

END 

***PRIM  P/^TH  DISTANCES’!'** 
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COMPUTER  PROGRAM  DEVERB 


* DEVERB  10/21/77  * 

* JEANIE  SAVAGE,  PROGRAMMER  * 

* SPECIFICATIONS  BY  RICK  BOSTIAN  =i‘ 

* IN  THIS  PROGRAM  DEVERBE RAT  I GN  * 

* IS  PERFORMED  IN  THE  FRE-  =*» 

* QUENCY  DOMAIN.  ♦ 

* LAST  CORRECTION:  12/06/77  * 


♦ 3it  4=  * 

INTEGER  XGRID 

INTEGER *2  I Z2 , 1 Z 3 , 1 Z A , I Z5 , 1 Z6 , YE S , I Z7 
DATA  YES/'Y  '/ 

DIMENSION  I STACK(20) , A( 1000) ,B( lOCO) , lARYI ICOO) 
DIMENSION  X(IOOO) ,Y(1000) 

DIMENSION  FINAL ( 1000 1 2) , IBEG( 20  ) , I END ( 2 0 ) , I CASE ( 20  ) 
COMMON  SF, I SF, IBM, THETA, SIGMA, GAMMA, C,D,RC0^FF,N2, 
&ISIZE,DB,OS,NBLK, IZ3,IZ4 

:fc  :Si  * :it  # :je  ^ * =(t  5{5  !(t # 5(5  >)t 

INITIALIZATION  ROUTINE 

* 5(1  :!t  :«c  * « ♦ ♦ 5)5  ♦ « 5(5  * 3(5  # ^(5  « 5): 

IREAD=0 
PRINT  500 

500  FORMAT! *0' , 'DEVERB'  ) 

PRINT  510 

510  FORMAT  ( 'O'  , 'NUMBER  OF  POINTS  PER  SIGNAL  (pCWER  OF  2)', 
• &•  (15)--') 

READ  520,  N2 

520  F0RMAT(I5) 

i.  PRINT  530 

530  FORMAT!'  ','IS  THIS  THE  FIRST  SIGNAL  OF  A MULTIPLE', 
£'  RUN?  !Y/N)-' ) 

READ  5A0, IZ2 

5A0  FORMAT !A1) 

PRINT  550 

550  FORMAT!*  ', 'SAMPLING  FREQUENCY  !F9.3) — ') 

READ  560,  SF 

560  F0RMATIF9.3) 

PRINT  570 

570  FORMAT!'  ', 'DIRECT  PATH  DISTANCE  !F9.5) — ') 

READ  580,  0 

580  F0RMAT(F9.5) 

PRINT  590 

590  FORMAT!'  ', 'SURFACE  PATH  DISTANCE  IN  METERS  (F9.5) — ') 
READ  580, DS 
PRINT  600 

600  FORMAT!'  'BOTTOM  PATH  DISTANCE  IN  METERS  (F9.5) — ') 
READ  580,  D8 
PRINT  610 

610  FORMAT!'  ', 'SURFACE  REFLECTION  TIME  IN  MSEC  (F9.5)  — 'j 
READ  580,  TS 
TS=TS/1000 
PRINT  620 

620  FORMAT!'  ', 'BOTTOM  REFLECTION  TIME  IN  MSEC  !F9.5) — ') 
READ  580,  TB 
TB=TB/1000 
PRINT  630 

630  FORMAT!'  ', 'BOTTOM  REFLECTION  CCEFFICIENT  (F9.5) — ') 
READ  58C,  RCOEFF 
WRITE!6,711) 

711  FORMAT!'  ','RMS  WAVE  HEIGHT  !F9.5) — ') 

READ  580,  SIGMA 
PRINT  720 

720  FORMAT!'  ', 'SURFACE  ANGLE  IF  INCIDENCE  UN  RADIANS)', 
&'  !F9,5)--') 

READ  580,  THETA 
PRINT  730 

730  FORMAT!'  ', 'BOTTOM  PHASE  SHIFT  (F9.5) — ') 

READ  580,  GAMMA 
PRINT  6A0 
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655 

660 

y 

\ 67  2 

67  5 

P 

\ . 

t, 

680 

690 

. 700 

i. 


1000 

1005 


lOlC 

102C 


1020 


1501 

1500 

1510 

1520 

1530 

1540 

1550 


FORMATC  *, 'SPEED  OF  SOUND  (F9.3)--') 

READ  560,  C 
PRINT  650 

FORMATC  'FREQUENCY  PLOT?  (Y/N)  — •) 

READ  540, IZ3 
PRINT  655 

FORMAT!'  't'TIME  PLOT  BEFORE  FFT?  (Y/N) — ') 

READ  540,177 

IF{ IZ7.NE.YES)  GOTO  659 

PRINT  657 

FORMATC  ', 'NUMBER  OF  POINTS  TC  BE  PLOTTED  (15) — ') 

READ  520,  INUMl 
PRINT  660 

FORMATC  ’,'TIME  PLOT  AFTER  FFT?  (Y/N)  — ') 

READ  540, IZ4 
1F( IZ4.NE.YES)  GOTO  672 
PRINT  657 
READ  520, INUM2 
PRINT  675 

FORMATC  ','ALL  BLOCKS?  (Y/N) — ') 

READ  540, IZ5 

IF(IZ5.EQ.YES)  GOTO  1000 

KPTR=0 

PRINT  680 

FORMATC  ', 'SPECIFIC  BLOCKS  (12)  (INPUT  99  WHEN', 

&'  FINISHED)—') 

READ  700,ITEMP 
FORMAT (I  2) 

IF( ITEMP.EQ.99)  GOTO  1000 
KPTR=KPTR+1 
ISTACK( KPTR)=ITEMP 
GOTO  690 

READ  DATA  TAPE 
^ * 

IF( IREAD.EQ.l ) GOTO  2000 
PRINT  1005 

FORMATC  ', 'READY  TO  READ  DATA  TAPE’) 

IF( 1Z2.NE.YES)  GOTO  1020 
READ  (5,1010)  KEND 
F0RMAT(8I6) 

DO  1030  1=1, N2, 8 

ITEMP=I+7  ; 

READ(5,1010)  ( lARY(J)  ,J  = I ,ITEMP)  i 

CONTINUE  I 

DETERMINE  DIVISION  BOUNDARIES  i 

♦ ’f' **=!=***  5i‘ =!= 

PRINT  1501 

FORMATC  '.'CONTINUING  WITH  CALCLLATIONS ' ) 

DO  1510  1=1, N2 
TIME=( I-l )/SF 
IF(TIME.LT.TS)  GOTO  1510 
ISF  = I 
GOTO  1520 
CONTINUE 
ISF=N2+1 
GOTO  1540 
DO  1530  1=1, N2 
TIME=( I-l )/SF 
IF(TIME.LT.TB)  GOTO  1530 
IBM=I 
GOTO  1550 
CONTINUE 
IBM=N2+1 

♦♦♦determine  NUMBER  OF  POINTS  IN  EACH  SECTOR*** 

NPT1=I SF-1 
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NPT2=  IBM-  ISF 

IF( IBM.LT.I SF)  NPT2=ISF-IBM 
NPT3=N2-1EM+1 

IF(  I8M.lt. ISF)  NPT3=N2-ISF+1 


C 

C 

C 

C 

C 

C 


C 

C 


C 

c 


c 

c 


c 

c 


* * !jr  * « * * ^«  # !(t  * 

DETERMINE  BLOCK  SIZE 

j|<  Jii  jjt  :j! 

***F1ND  SMALLEST  AND  MIDDLE  NUMBER  OF  POIN 
1700  ISMALL=MINC(NPT1 ,NPT2) 

ISMALL  = MIN0(NPT3, I SMALL) 

IF( ISMALL.EC.NPTl ) MI DDLE=M INO( NPT 2t NPT 3 ) 
IF(  ISMALL.EC.NPT2  ) M I DDLE=M  INO ( NPTl , NPT 3 ) 
IF( ISMALL .EC.NPT3 ) MI CDLE=M  INO ( N PT 1 ,NPT 2 ) 

***DET'feRMINE  NUMBER  OF  POINTS  PER  BLOCK*** 
2000  IF( lABSI Ml DDLE/2-TSMALL).GT. (MIDDLE- I SMALL 
CISIZE=ISMALL 

IF( ISIZE.EQ.ISMALL)  GOTO  2200 
DO  2010  K=l,30 

IF( ISMALL.EC.MIDDLE/K ) GOTO  2020 
IF( ( ISMALL.EQ.MIDDLE/K) .LT. ( I SMAL L-MI DDL E/ 
CGCTO  2020 
2010  CONTINUE 
2020  ISIZE=MIDOLE/K 

***DETERMINE  NUMBER  OF  BLOCKS  PER  SECTOR** 
2200  N8LK1=NPT1/ISIZE 
N8LK2=NPT2/ISI ZE 
NBLK3=NPT3/ISIZE 

***DETERMINE  NUMBER  OF  POINTS  SKIPPED  EACH 
2300  NSKI P1=NPT1-NBLK1*I SIZE 
N SKI P2=NPT2-NBLK2*I SIZE 
NSKIP3=NPT3-NBLK3*ISI ZE 

***PRINT  SECTOR  INFORMATION*** 
WRITE(6,2A05)  I SI ZE 
2405  FORMAT ( 'O' f • ISIZE=' , 13) 

IT1=1 
IT2  = 2 
IT3  = 3 

WRITE(6 t241C)  I T1 tNBLKl ,NSK IPl 
2410  FORMAT! 'O' , 'SECTOR  'til,'  CONTAINS  •tI2,' 
$2X,I3,'  POINTS  SKIPPED') 

WRITE(6»241C)  I T2  , NB LK2 » NSK I P2 
WRITE(6,241C)  I T3 » NBLK3 , NSK I P3 


C 

C BUILD  STACK  IF  PROCESSING  ALL  BLOCKS 

c 

2500  IF( IZ5.NE.YES)  GOTO  2700 
NUMBLK=N2/I SIZE 
DO  2510  1=1 tNUMBLK 
ISTACK( 1 )=I 
2510  CONTINUE 

KPTR=NUMBLK 

C 

C DETERMINE  BEG  6 END  POINTS  £ CASE 

C 

2700  16=0 

C 

C ♦**BLOCKS  IN  1ST  SECTOR*** 

ITEMP=NSKIP1+1 
DC  2710  I=1,M3LK1 
16=16+1 

IBEG( 16  )=ITEMP 

IEND( I6)=ITEMP+ISIZE-1 
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TS  IN  SECTOR 

) ) 

( K+1)  ) ) 

* 

SECTOR*** 


BLOCKS, ' , 


ICASE( I6)=l 
ITEMP=I END( I6)+l 
271C  CONTINUE 
C 

C *=»*BLOCKS  IN  2ND  SECTOR*** 

ITEMP=I SF 
ITEMP1=NSKI P2 
DO  2720  I*=1,NBLK2 
16=16+1 

IBEG(  16  )=ITEMP 
IEND( I6)=ITEMP+ISIZE-1 
ITEMP=I END(  I6)  + l 
ICASE( I6)=2 

IF(  ISF.GT.IEMJ  ICASE(T6)  = 3 
IF( ITEMP2-EG.0)  GOTO  2720 
ITEMP=I TEMP+1 
ITEMP1=ITEMP1-1 
272C  CONTINUE 
C 

C ***BLOCKS  IK  3RD  SECTOR*** 

DO  2730  I=1,NBLK3 
16=16+1 

IBEG(  I6»=IBK+(  I-l  )*1SIZE 
IF(  ISF.GT  .1 BM)  IBEG( I6)  = ISF+( I-1)*ISIZE 
IEND(  16  )=1BEG(  16  )+ISI ZE-1 
ICASE( I6)=4 
273C  CONTINUE 
C 

C PRINT  TIME  PLOT  BEFORE  FFT 

C 

290C  IF( IZ7.NE.YES)  GOTO  3000 

IF{ INUMl.GT.lOOO)  GOTO  2910 
NUM  = I NUMl 
KPL0TS=1 
GOTO  2920 

2910  KPL0TS= INUMl/1000+1 
NUM=1000 

2920  DC  2940  I=1,KPL0TS 

IF( { I.EQ.KPLOTS) . AND. (I .NE. 1) ) NUM=I NUM 1 - ( KFLCTS- 1 ) 
&*1000 

DC  2930  J=1,NUM 

Y( J)=FLOAT( lARYl (I-l) *1000+J) ) 

X( J>=( ( I-1)*1000+J-1)/SF 
2930  CONTINUE 

XGRI0=NUM/1C 

IF(MOD( NUM, 10).ME.0)  XGRI D=XGR I C+ 1 
CALL  PL0TDV(X,Y,NUM,XGRI0,3 ,NBLK) 

2940  CONTINUE 
C 

C PROCESS  BLOCKS 

c 

3000  DO  3110  1=1  ,KPTP 
NBLK=ISTACK(n 
C 

C ***DETERMINE  BOUNDARIES  AND  CASE  CF  BLOCK*** 

I BEGO=I eEG( N6LK ) 

IENDO=I END( KBLK ) 

ICASEO=  ICASE(NBLK  ) 

C 

C ***PROCESS  BLOCK*** 

PRINT  3015,  IBEGO,IENDO,ICASEO 
3015  FORMAT!///'  • , • I B EG= • , 14, 4X , * I END= * , 14 , 4X, ' IC ASE= ' , 
K = 1 

DO  3020  J=I BEGO, I ENDO 
A(K)  = FLOAT(  IARY( J ) ) 

B(KJ=0.0 
K = K+1 


ID 
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C/SLL  SIGNAL  ( A,  B , F INAL  , I BEGO  , I C A S EO  , I S I Z E ) 

IFKI  .NE.l)  .OR.(NSKIPl.EQ.O))  GGTC  30A0 
DC  3030  J=1,NSKIP1 

SL0PE=FINAL(IBEG(!)+1,1)-FINAL( IBEG(I),1  ) 

F INAL( J,1 )=-SLOPE+FINAL( I8EG( I ) , 1 ) 

SL0PE  = FINAL( IBEG( I ) +I , Z ) -F I N AL ( I B EG ( I ) , 3 ) 

F INAL(  J»2)=-SL0PE+F1NAL(IBEG(  1)  ,2) 

3030  CONTINUE 

3040  IF( ( I .GT.NBLK1+NBLK2 ) .OR. ( I .LE.NELKl+l ) ) GOTO  3110 
IFINSKI P2.EC.0)  GOTO  3110 
DO  3050  J=1,NSKIP2 

ite:^p=i  end(  nblkh-i  ) + i 

TEMP=FI NAL ( ITEMP+1 ,1 )-FINAL ( ITEMP-1, 1 ) 
FINAL(ITEMPa)  = FINAL(  I TEMP- 1 , 1 ) +T  EMP 
TEMP=FINAL( ITEMP  + 1 , 2 ) -F I N AL ( I T EM  P-1 , 2 ) 

FINAL! I TEMP, 2)  = F INAL ( ITEMP- 1 , 2 ) +T EMP 
3050  CONTINUE 
3110  CONTINUE 
C 

C PRINT  TIME  PLOT  AFTER  FFT 

C 

3300  IF( IZ4.NE.YES)  GOTO  3500 
DO  3308  I=1,KPTR 
NBLK=ISTACK{I ) 

N=ISIZE 

I BEG0=I BEG( NBLK ) 

IEN00=I END(NBLK) 

K=1 

DC  3304  J=IBEG0,I FNDO 
A (K)=F INAL! J,1 ) 

B(K)=FI NAL! J,2) 

K=K+1 

3304  CONTINUE 

CALL  CPFT2! A,B,N ,N,N, 1 ) 

DO  3306  J=1  ,ISI ZE 
FINAL!IBEG0+J-1,1)=A! J) 

3306  CONTINUE 
3308  CONTINUE 

IFINSKIPl .EC.O)  GOTO  3313 
DC  3312  I=1,NSKIP1 

FINAL!  I ,l)  = !FINAL!ISEG!l),l)/(NSKIPl  + l))=i-I 

3312  CONTINUE 

3313  IF!NSKIP2.EC.0)  GOTO  3315 
DO  3314  I=1,NSKIP2 
ITFMP=! NBLKl+I )+l 

TEMP=FI  NAL!  ITEMP<-1 , 1 )-FINAL  ! ITEMP-1, 1 ) 
FINAL(ITEMP,1)=FINAL! I TEMP-1 , 1 ) +T EMP 

3314  CONTINUE 

3315  IF! INUM2.GT.IEND! KPTR  ) ) I NUM2= I END  ! KPTR ) 
IF!INUM2.GT.1000)  GOTO  3310 
NUM=INUM2 

KPL0TS=1 
GOTO  3320 

3310  KPL0TS=INUM2/1000+1 

NUM=1000 

3320  DO  3340  I=ltKPLOTS 

IF! ! I .EQ.KPLOTS ) .AND.  ! I .NE. 1) ) NUM  = INUM2- ! KFLCTS-1) 
8*1000 

DC  3330  J=1,NUM 
Y! J)=FINAL( !I-1 ) *1000+J, 1 ) 

Y!J)  = Y! J) /FLOAT ! I SIZE) 

X!J)=! ! I-l )*1000+J-1)/SF 
3330  CONTINUE 

XGRID=NUM/1C 

IF! MOO! NUM,10).NE.0)  XGRI D=XGR I 0 + 1 
CALL  PL0T0V!X, Y,NUM,XGRID,2 ,NBLK  ) 

3340  CONTINUE 

r 

C CONCLUSION 
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C 

3500  PRINT  3510 

3510  FCRMAK  • ','ARE  YOU  FINISHED?  (Y/N)—') 

READ  540, I Z6 
IFnZ6.  EQ.YES)  STOP 
IREA0=1 
GOTO  672 
UE3UG  SUBCHK 
tND 
C 

c 

C ****SIGNAL  SUBRCUTINE*=!=’!'>:' 

C 

SLBROUT  INE  S IGN AL ( A , B , F I NAL , IBEGO  , ICASEO , N ) 
DIMENSION  FINAL(1000,2) 

DIMENSION  X{ 1000) , Y( 1000) ,A( 1000)  ,B( 1000) 

INTEGER*?  YES,IZ3,IZ4 

COMMON  SF,  ISF,  I BM,THFTA,SIGMA,GAMNa,C,D,RC0E'='",N2, 
{1ISIZE,D8,0S,NBLK,IZ3,IZ4 
INTEGER  XGRID 
REAL  KO 

DATA  YES/*Y  •/ 

C 

C PERFORM  FFT 

Q *X«5!<******** 

C 

100  CALL  CFFT2(A,B,N,N,N, -1) 

C 

C PERFORM  CORRECTIONS 

c 

500  DO  520  1=1, N 

FINAL! I 8EG0+I-1 , 1 )=A( I ) 

FINAL! I 3EG0  + I-1 ,2  ) = B(  I ) 

FREQ=!I-1)*SF/N 
PI=3. 14159 
K0=2.*P I*FREQ/C 
IF! ICASEO. EQ.l ) GOTO  520 
IF!  ICASEO. EG. 3)  GOTO  510 
C 

C ***CORRECTICN  FOR  SURFACE  REFLECTION*** 

G=! !4.*PI*SIGMA*FPE0/C)*C0S!THETA)  )**2 
ITEMP=(  IBEGC+I-1 ) -!  I SF-1 ) 

TMAG= SORT (FINAL! I TEMP , 1 ) **2 +FI N AL  ! ITEMP , 2)**2  ) 
TPHASe=ATAM2( FINAL! ITEMP, 2 ) , F I N A L ! IT EMP , 1 ) ) 
S=D/DS*EXP! -G/2. ) *TMAG 
SPHASE=TPHAS£-! DS-D)*KO-PI 

FINAL! IBEGO  + I-1 , 1 )=FINAL! I B EGO+ I - 1 , 1 ) -S *COS (SPHASE  ) 
FINAL! I BEGO  + I-1 , 2 ) = FINAL!  IB EGO  + I - 1 , 2 ) -S*S I N ! S PH AS E ) 
IF! ICASEO. EG. 2)  GOTO  520 
C 

C ***CORRECTION  FOR  BOTTOM  REFLECTION*** 

510  ITEMP=! IBEGC+I-1 )-!IBM-l) 

TMAG=SQRT (^INAL ! I TEMP , 1 ) **2 + FIN AL ! IT EWP , 2 )**2  ) 
TPHASE=ATAN2! FINAL! ITEMP, 2) , FINAL! ITEMP , 1 ) ) 
S=RCOEFF*D/CB*TMAG 
SPHASE=TPHASE-! DB-D)*K0+GAMMA 

FINAL! IBEGO  + I-1, 1 )=FIMAL( 1 8 EGO+ 1 - 1 , 1 ) -S*C CS ( SPHAS E ) 
FINAL! IBEGO  + I-1 ,2 )=FINAL! IB EGO+I- 1,2 )-S*S  IN  (SPHASE) 
520  CONTINUE 
C 

C PRINT  FREQUENCY  AND  TIME  PLOTS 

C 

1000  IF! IZ3.NE.YES)  GOTO  1020 
C 

C ***FREQUENC r PLOT*** 

K = 0 
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ITEMP=N/2 
DO  1010  I=l,ITEMp 
K=K  + 1 

Y(K)=FINAL ( IBEG0+1-1» l)*-2+FTNAL ( I BEGO+ I- 1 , 2 ) 
Y(K)=10*ALOG10( Y(K>*0) 

X(K)=(K-1)*SF/N 
1010  CONTINUE 
NUM  = K 

XGRI0=NUM/1  C 

1F(M00( NUM, 10) .NE  .0)  XGRIO=XGRI C + 1 
CALL  PL0T0V(X, Y,NUM,XGRID, 1 ,N8LK) 

1020  RETURN 

DEBUG  SUBCHK 
END 
C 
C 

c 

Q SUBROUTINE’!'*** 

C 

SUBROUTINE  PLOTOV ( X» Y , N, XGR I 0, M , NE  ) 

INTEGER  OtXCRID, YGRIOtAXIS 

D I MENS  I ON  Y(1000),C7(101),0(6)tX(  1000) ,KAXI  S( 51) 
DATA  lOASH/lH-/ , I STAR /I H*/,  lOOT/lH./ 

DATA  IBAR/lHI/^ IPLUS/1H+/, I BLANK/lh  /,IX/1HX/ 
AXIS=51 
YGRID=6 
XGRID=XGRIO-H 
2120  N1=N-1 

Y6=Y(1 ) 

Y1=Y( 1 ) 

DO  2200  1 = 1 , N 

I F(Y6-Y( I ) .GE.0.0  ) GOTO  2180 
Y6  = Y( I ) 

2180  1F(Y1-Y( I ).LE.O.O)  GOTO  2200 

Y1=Y( I ) 

2200  CONTINUE 

S=Y1*(AXIS-1)/(Y6-Y1 ) 

X1=X(1) 

Xia=X(N) 

0( 1)=Y1 

0(6)=Y6 

IIX=XGRID-1 

DO  2410  1 = 1, IIX 

C7( 1)=X((I-1)*10+1) 

2410  CONTINUE 

C7(XGRI D)=C7(XGRID-1 ) ^10* ( X ( 2 )-X ( 1 ) ) 
IF(N.EQ.(XGRI0-1 )*10)  C7(XGRID)=X(N) 

I IY=YGRID-1 
DC  2440  1 = 2, IIY 

0(1 ) = (FLOAT( I-l )* (Y6-Y1 ) / FLO  AT ( YGR  I D- 1 ) ) +Y1 
2440  CONTINUE 

WRITE( 6,2460) 

2460  FORMAT!///,*  •) 

IF(M.NE.l)  GOTO  2485 
PRINT  2470, NB 

2470  FORMAT! • 0 ' , 32X, * 3 LOCK • , IX , I 2 ) 

2485  IF(M.EQ.l)  PRINT  2486 

2486  FCRMATI*  •,27X,'DB"S  VS.  FREQUENCY') 

IF(M.EQ.2)  PRINT  2488 

2488  FORMATC  ',  23X,  • VOLT  AGE  VS.  TIME  (AFTER  FFT)') 
IF(M.EQ.3)  print  2487 

2487  FORMATC  *,  23X ,' VOLTAGE  VS.  TIME  (BEFORE  FFT)') 
WRITE(6,2500)  ( 0(  1) , I = 1 , YGR I D ) 

2500  F0RMAT(9X,11(1PE10.2)  ) 

S1=(X!0-X1)/10.0*(XGRID-1) 

D = 1 
Ll=l 
L = 1 

IZ=IFIX(-S+1.5) 

ITEMP=(  XGRIO-1  )*10'-1 
DO  2900  I1=1,ITEMP 
IF(N.LT.Il)  GOTO  2510 
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YTEMP=(  Y(  !1)*FL0AT(YGRID-1)  ’t'lO.  0 / (Y6-Y1 ) )-S 
I Y=IFIX(YTFVP+1 .5) 

2510  IFILI.GT.LJ  GOTO  2760 
00  2650  IP=1,AX!S 
KAXIS( IP)=IDASH 
2650  CONTINUE 

DO  2680  1 = 1 , AXIS,  10 
KAXIS( I )=IPLUS 
2680  CONTINUE 

IF(N.LT.Il)  GOTO  2720 

IFKYl.LE.O.O)  .AND.(0.0.lE.Y6))  K AXI  S ( I Z ) = I COT 
KAXIS( IY)=I STAR 

2720  WRITE(6,2725)  C 7 ( D ) , ( KAX I S ( J ) , J = 1 , AX ! S ) 

2725  FORMAT! 1PE13. 2, 2X,115A1) 

L1=L1+10 
D = 0+1 
GCTO  2870 

2760  DO  2780  IP=1,AXIS 
KAX1S( IP)=1ELANK 
2780  CONTINUE 

DC  2810  1 = 1, AXIS, 10 
KAXISI I )=IBAR 
2810  CONTINUE 

IF(N.LT,I1)  GOTO  2860 

IF((Y1.LE.0.0).AND.(0.0.LE.Y6))  KAXISdZ  ) = ICGT 
KAXTS( IY)=ISTAR 

2860  WRITE(6  ,2  865)  ( KAXIS ( J ) , J = 1 , AX  I S ) 

2865  FORMAT!  15X,  115A1  ) 

287C  L=L+1 

IF!L.GT.!XGRID-l)*10+2)  GOTO  2910 
2900  CONTINUE 
2910  RETURN 

CE8UG  SU8CHK 
END 
C 
C 

C ***FOURIER  TRANSFORM  SUBROUT  INE=^=«=* 

C 

SUBROUTINE  CFFT 2 ! A , B , NTOT , N , NSP AN  , I SN ) 


I 
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COMPUTER  PROGRAM  TDEVERB 


50C 

510 

520 

540 

550 

560 

570 

580 

590 

600 

610 

620 

630 

711 

720 

730 

640 

645 


* TDEVERB  12/07/77  * 

=*  JEANIE  SAVAGE,  PRCGRAMMER  * 

* SPECIFICATIONS  BY  RICK  BOSTIAN  « 

« IN  THIS  PROGRAM  DEVER6ERATICN  * 

* IS  PERFORMED  IN  THE  TIME  * 

* DOMAIN.  * 

* LAST  CORRECTION:  12/08/77  =!= 


INTEGER  XGRID 

INTEGER*2  I 2 2,  I Z 3 , 1 Z 4 , ! Z5, 1 Z6, Y E S , I Z 7 
DATA  YES/'Y  •/ 

DIMENSION  I STACK (20) , A ( 1000 ) , S ( 1 000 ) , I ARY ( 1 COO ) 
DIMENSION  X(IOOO) , Y( 1000) 

DIMENSION  FINAL! 1000) ,IBEG(  20 ) , I END ( 20 ) , I C ASE ( 20 ) 

:!>;!< *#«*  3|:« 

INITIALIZATION  ROUTINE 


IREAD=0 
PRINT  500 

FORMAT { 'O' , 'DEVERB* ) 

PRINT  510 

FORMAT( *0* , 'NUMBER  OF  POINTS  PER  SIGNAL  (PCWER  0=  2)y 
C (15)  — ') 

READ  520,  N2 

FaRMAT(I5) 

FORMAT (A1 ) 

PRINT  550 

FORMATC  *, 'SAMPLING  FREQUENCY  (F9,3) — ') 

READ  560,  SF 

FQRMAT(F9.3) 

PRINT  570 

FORMATC  ', 'DIRECT  PATH  DISTANCE  (F9.5)  — ') 

READ  580,  0 

F0RMAT(F9.5) 

PRINT  590  . 

FORMATC  ', 'SURFACE  PATH  DISTANCE  IN  METERS  (F9.5)~-'J 
READ  580, DS 
PRINT  600 

FORMATC  ', 'BOTTOM  PATH  DISTANCE  IN  METERS  (F9.5)  — ') 
READ  580,  DB 
PRINT  610 

FORMATC  ', 'SURFACE  REFLECTION  TIME  IN  MSEC  (F9.5)  — •} 
READ  580,  TS 
TS=TS/1 000 
PRINT  620 

FORMATC  ', 'BOTTOM  REFLECTION  TIME  IN  MSEC  (F9.5) — ') 
READ  580,  TB 
TB=TB/1000 
PRINT  630 

FORMAT!'  ', 'BOTTOM  REFLECTION  COEFFICIENT  (F9.5) — •) 
READ  580,  RCOEFF 
WRITE!6,711  ) 

FORMAT!'  ','RMS  WAVE  HEIGHT  !F9.5)  — ') 

READ  580,  SIGMA 
PRINT  720 

FORMATC  ', 'SURFACE  ANGLE  OF  INCIDENCE  (IN  RADIANS)', 
£'  (F9.5)--') 

READ  580,  THETA 
PRINT  730 

FORMAT!'  ', 'BOTTOM  PHASE  SHIFT  (F9.5) — •) 

READ  580,  GAMMA 
PRINT  640 

FORMAT!'  ', 'SPEED  OF  SOUND  !F9.3)~') 

READ  560,  C 
PRINT  645 

FORMAT!'  ', 'BLOCK  SIZE  (15)  — ') 

READ  520,  ISIZE 
PRINT  650 
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650  FORMATC  'FREQUENCY  PLOT  BEFORE  CORR  ECT I CNS?  • , 

£•  (Y/M) — ') 

REAO  540, IZ3 
PRINT  652 

652  FCRMATI • 'FREQUENCY  PLOT  AFTER  CORRECTIONS?  lY/N) — •) 
READ  540,  IZ2 
PRINT  655 

655  FORMAT!'  ','TIME  PLOT  BEFORE  CORRECTIONS?  (Y/MJ — ') 
READ  540, IZ7 
IF! IZT.NE.YES)  GOTO  659 
PRINT  657 

657  FORMAT!'  ', 'NUMBER  OF  POINTS  TO  BE  PLOTTED  115) — •) 
READ  520,  INUMl 

659  PRINT  660 

660  FORMAT!'  ','TIME  PLOT  AFTER  CORRECTIONS?  (Y/N) — •) 
READ  540, IZ4 

IF!  IZ4.NE.YES)  GOTO  672 
PRINT  657 
READ  520,INUM2 
672  PRINT  675 

675  FORMAT!'  •,»ALL  BLOCKS?  (Y/N) — •) 

READ  540, IZ5 

IF! IZ5.EQ.YES)  GOTO  1000 

KPTR=0 

PRINT  680 

680  FORMAT!'  ', 'SPECIFIC  BLOCKS  !I2)  (INPUT  99  WFEN  FIMS 
690  READ  700,ITEMP 
700  FCRMAT(I2) 

IF(  1TEMP.EQ.99  ) GOTO  1000 
KPTR=KPTR+1 
ISTACK! KPTR  )=ITEMP 
GOTO  690 

:(C :{t  * ijt  * it  * 

READ  DATA  TAPE 

it  irt)t«it!j<ititit5fti:iti:is 

1000  IF! IREAO.EQ.l ) GOTO  2500 
PRINT  1005 

1005  FORMAT!'  ', 'READY  TO  READ  DATA  TAPE') 

1010  F0RMAT(8I6) 

1020  DO  1030  1 = 1, N2, 8 
ITEMP=I+7 

READ(5,1010)  ( lARY! J)  ,J  = I ,ITEMP) 

1030  CONTINUE 

it  it  it  it  it  1}:  * it  it  it  it  tv  it  A it  it  * it  it  it  it  it  it  * * 

DETERMINE  DIVISION  BOUNDARIES 

it  it  it  it  it  it  it  Ait  it  A it5)t  it  i it  it « it  it  a a it  it  it  it  it  i it 

PRINT  1501 

1501  FORMATC  ' , ' CONT  INUl  NG  WITH  CALCULATIONS') 

1500  DO  1510  1=1, N2 
TIME=(I-1 )/SF 
1F1TIME.lt. TS)  GOTO  1510 
ISF=I 
GOTO  1520 
1510  CONTINUE 
ISF=M2+1 
GOTO  1540 

1520  CO  1530  1=1, N2 
TIME=(I-1)/SF 
IFITIME.LT.TB)  GOTO  1530 
IBM=I 
GOTO  2500 
1530  CONTINUE 
1540  IBM=N2+1 


itititititititititititititititititititititititititititiiititititirit 

BUILD  STACK  IF  PROCESSING  ALL  BLOCKS 

it  ♦ it  * it  it  it  it  * it  it  it  * it  it  it  it  it  it  it  # it  it  it  it  it  * it  it  it  it  it  it  it  it 


c 

2500 


251C 

C 

C 

C 

c 

c 

2900 


IF(  IZ5.NE.YFS)  GOTO  2900 

NUM8LK=N2/I SIZE 

DO  2510  I=1,NUMBLK 

ISTACKI I )=I 

CONTINUE 

KPTR=NUM8LK 

PRINT  TIME  R.OT  BEFORE  CORRECTIONS 

:(c  Jji * :(c :j(  :(!##* 5)t:{s  :(< j^s  jjt :)[ 3{t 3j!  :jt  jjt ^ ^ s(i * :jc :jt ^ 3)t 


IF(  IZ7.NE.YES) 
IF( IZ5.NE.YES) 


GOTO  3000 
GOTO  2950 


2910 

2920 


2930 


2940 


2950 


296C 


2970 

C 


♦ ♦♦time  plot  fop  ENTIRE  SIGNAL=i=** 

IF( INUMl.GT.lOOO)  GOTO  2910 

NUM=INUM1 

KPLOTS= 1 

GOTO  2920 

KPLOTS=INUM1/1000+1 

NUM=1000 

DC  2940  1=1  ,KPLOTS 

IF( < I .EQ.KPLOTS ) . AND. {I .NE. 1) ) NUM=INUMl-( KFLCTS-1 )* 
CICOO 

DC  2930  J=1,NUM 

Y( J»=FLOAT( 1ARY( ( 1-1 ) *1000+ J J ) 

X(J)=( ( I-l)*1000+J-l) /SF 

CONTINUE 

XGR1D=NUM/1 0 

IFlMODt  NUM, 10) .ME .0)  XGR10=XGRI C+1 
NBLK=-I 

CALL  PL0T0V(X,Y,NUM,XGRID,3,NBLK) 

CONTINUE 
GCTO  3000 

***TIME  PLOT  FOR  INDIVIDUAL  BLOCKS*** 

) DO  2970  I=1,KPTR 
NELK=ISTACK (I ) 

DO  2960  J=1,ISIZE 

Y(J)=FLOAT( IARY(NeLK-l)*ISI ZE+J) 

X( J)=( ( NBLK-1 )*ISI ZE+J-1 )/SF 
3 CONTINUE 
XGRID=ISIZE/10 

IF(MOD(  IS  IZE,10) .NE.O)  XGRI D=XGR I C + 1 
CALL  PLCTDVIX, Y, I SIZE,XGRI0,3,NBLK) 

CONTINUE 


♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

PERFORM  CORRECTIONS 
♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦* 

'3000  DC  3010  1 = 1, N2 

FREQ= (I-l )*SF/ISI ZE 
PI=3. 14159 

G=( (4.*PI*SIGMA*FREQ/C)*C0S(THETA)  )**2 
FINAL!  I I=FLCAT(1ARY(  I)  ) 

ITEMP=I-ISF+1 

ITEMP2=I-I3M+1 

IF(I.GE.ISF)  FINAL! I ) =F IN AL ! I ) + EX F ! -G/ 2 . ) *D*DS* 
&FINAL!ITEMP) 

IF!  I. GE  .IBM)  FINAL (I )=FINAL ! I )-RCCEFF*D*DB* 
SFINAL! ITEMP2) 

.3010  CONTINUE 

♦♦♦♦♦A ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦ 

PRINT  TIME  PLOT  AFTER  CORRECTIONS 
♦♦♦♦♦♦♦♦♦*♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

IF!  IZ4.NE  .YES)  GOTO  3400 
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IF(U5.NE.YES)  GOTO  3350 
C 

c *«4.TIME  plot  for  entire  SIGNAL=<‘*=«' 

IF(  INUM2.GT.100C)  GOTO  3310 

NUM=INUM2 

KPL0TS=1 

GOTO  3320 

3310  KFL0TS=INUM2/1000+1 
NUM=1000 

3320  OG  3340  I=1,KPL0TS 

IF( ( I .EQ.KPLOTS ) .AND. (I .NE. 1) ) NUW=INUM 2- < KPLCTS- 1 ) * 
&1COO 

DC  3330  J=1,NUM 

Y{ J>=FI  NAL( ( I-l )*1000+J) 

X(J)  = ((  I- 1) *1000  + 0-1 )/SF 
3330  CONTINUE 

XGRI0=NUM/1 G 

IF(MOO( NUM, lOJ.NE.O)  XGRIO=XGRI C+1 
N8LK=-1 

CALL  PL0TDV(X,Y,NUM,XGRIDf2 ,NBLK) 

3340  CONTINUE 
GCTQ  3400 
C 

C ***TIME  PLOT  FOR  INDIVIDUAL  BLOCKS*** 

3350  DO  3370  I=1,KPTR 
NeLK=ISTACK(n 
00  3360  J=1,ISIZE 
Y(JI= FINAL! <NBLK-1)*ISIZE+J) 

Xt J)=( { N3LK-1)*ISIZE+J-1)/SF 
3360  CONTINUE 

XGR1D=1SIZE/10 

IF(MOD( ISIZE.IO) .NE.O)  XGRI D=XGR I C+1 
CALL  PLCT0V(X,Y,ISIZE,XGRID,2tNELK) 

3370  CONTINUE 
C 
C 

C PRINT  FREQUENCY  PLOTS 

C 

34C0  IF( IZ3. NE.YES)  GOTO  3450 
C 

C ***FREQUENCY  PLOT  BEFORE  CORRECTIONS*** 

DO  3430  1=1  ,KPTR 
NBLK=ISTACK(n 
DC  3410  J=1,ISIZE 

A( J)=FLOAT( IARY( ( N8LK-1 » * IS  I Ze+ J > ) 

B( J)=0.0 
341C  CONTINUE 

CALL  CFFT2( A,8,IS  IZE, ISIZE, ISI ZE ,-l) 

ITEMP=I SI ZE/2 
DO  3420  J=1jITFMP 
Y(  J)=A( J)**2+S ( J )**2 
Y( J)  = 10*AL0G10(Y( J)  ) 

X( J)=( J-1 )*SF/ISIZE 
3420  CONTINUE 

XGRID=ITEMP/10 

IF  (M00(  ITEMP, 10 ) .NE.O)  XGR ID=X GR I C+1 
CALL  PL OTOV(X,Y, ITEMP  ,XGRI0 , 1 ,N ELK ) 

3430  CONTINUE 
C 

C ***FREQUENCY  PLOT  AFTER  CORRECTIONS*** 

3450  IF(IZ2. NE.YES)  GOTO  3500 
DO  3480  I=1,KPTR 
NeLK=ISTACK(I ) 

DC  3460  J=1,ISIZE 

AIJ)  = FINAL1 (NBLK-1 )*ISIZE+J  ) 

B(  J)  = 0.0 
346C  CONTINUE 

CALL  CFFT2( A,Bt ISIZE, ISIZE, ISI 2E,-1) 

ITEMP=I SIZE/2 
DO  3470  J=1,ITEMP 
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Y(  J)=A(  J)=!‘*2+B(  J)**7 
Y(  J)  = 10.=!‘ALCG10(  Y(  J)  ) 

X(  J)=(J-U*SF/ISIZF 
347C  CCNTINUE 

XGRID=I TEMP/10 

IF(MQD(ITEMP,10).NE.0)  XGRI D=XGR I C+1 
CALL  PLOTOV(X,Y,!TEMP  ,XGR I D , 4 , N E LK ) 

348C  CCNTINUE 
C 

C CONCLUSION 

C jji#*** 

c 

3500  PRINT  3510 

3510  FORMAT! • '.'ARE  YOU  FINISHED?  (Y/N)--') 

READ  540,  IZ6 
IF! IZ6. EQ.YES)  STOP 
IREAD=1 
GCTO  672 
DEBUG  SUBCHK 
END 
C 
C 

c 

c 

C SU8ROUTINE=4=*^* 

C 

SUBROUTINE  PLOTDV ! X , Y , N, XGR I D, M , NE  ) 

INTEGER  D,XGRID,YGRIO,AXI S 

DIMENSION  Y!1000) , C7 ! 101 ) , 0 ! 6 ) , X ! 1000) ,KAXI S! 51) 

DATA  IDASH/IH-/, I STAR/IH*/, I DOT /IF./ 

DATA  I3AR/1HI/,IPLUS/1H+/,IBLANK/1H  /,IX/1H>(/ 

AXIS=51 
YGRID=6 
XGRI0=XGRID+1 
2120  N1=N-1 

Y6=Y!1) 

Y1=Y!1) 

00  2200  1=1  , N 

IF!Y6-Y! I ) .GE.0.0)  GOTO  2180 
Y6=Y! I ) 

2180  I F! Yl-Y! I ) .LE.0.0  ) GOTO  2200 
Y1=Y! I ) 

2200  CONTINUE 

S = Y1*!  AXIS-1)/!Y6~Y1) 

X1=XI 1 ) 

X10=X(N) 

0! 1)=Y1 

G!6)=Y6 

I IX=XGRID-1 

DC  2410  1=1, IIX 

C7!I)=X!!I-1)*10+1) 

2410  CCNTINUE 

C7!XGRID)=C7!XGRID-1 ) +10* ! X ! 2) -X  ! 1 ) ) 

IF1N.EQ.!XGRI0-1 )*10)  C7! XGRIO) =X( N) 

I IY=YGRID-1 
DC  2440  1=2, IIY 

0!  I )=!  FLO  ATll-1  )*!Y6- YD/FLO  AT!YGR  ID-1)  )+Yl 
2440  CONTINUE 

WRITE! 6,2460) 

2460  FORMAT!///,'  ') 

IF!NB.GT.O)  GOTO  2466 
Ne=-NB 

PRINT  2465,  N8 

246  5 FORMAT! ' 0 • , 32X, • PLOT • , IX, 1 2 ) 

GOTO  2485 

2466  PRINT  2470, NB 
2470  FORMAT! 'O' , 32X , ' B LOC K ’ , IX , I 2 ) 

2485  IF!M.EQ.l)  PRINT  2486 

2486  FORMAT!'  •,17X,'DB"S  VS.  FREQUENCY  !BEF0RP  CCRRECTICN 
IF(M.EQ.2)  PRINT  2488 

2438  FORMAT!'  ', 20X, * VOLTA GE  VS.  TIME  !AFTER  CORRECTIONS)') 


46 


IF(M.EQ.3)  PRINT  2487 

2487  FQRMATC  20X  ,•  VOLTA  GE  VS.  TIME  (BEFORE  COPR  ECT  I CNS  ) ') 

IF(M.EQ.4)  PRINT  2489 

2489  FCRMATC  S17X,'Oe"S  VS.  FREQUENCY  (AFTER  •» 
CCORRECTIONS)'  ) 

WRITE(6,2500)  ( 0 ( I ) , I =1 » YGR  I D ) 

2500  F0RMAT(9X,11(1PE10.2)  ) 

S1=(X10-X1 )/10.0*(XGRID-l ) 

0=1 
Ll  = l 
L = 1 

IZ=IFIX(-S+1.5) 

ITEMP=( XGRIC-1 )*1C+1 
DO  2900  Il=l,ITEMP 
IF(N.LT.Il)  GOTO  2510 

YTEMP=(Y(  II  )*FL0AT(YGRI0-1) ’«'10.0/(Y6-Y1 ) )-S 
IY=IFIX(YTEMP+1.5) 

2510  IF(Ll.GT.L)  GOTO  2760 
00  2650  IP=1,AXIS 
KAXIS( IP)=IDASH 
2650  CONTINUE 

DC  2680  1=1, AXIS, 10 
KAXIS( I )=IPLUS 
2680  CONTINUE 

IF(N.LT.Il)  GOTO  2720 

IF( (Y1.LE.0.0).AN0.(0.0.LE.Y6) ) KAXI S ( I Z ) =I COT 
KAXIS( IY)=I STAR 

2720  WRITE(6,2725J  C 7 ( 0 ) , ( KAX I S ( J ) , J = 1 , AX  I S ) 

2725  FORMAT! 1PE12. 2, 2X,115A1J 
L1=L1+10 
0 = 0+1 
GCTO  2370 

276C  DO  2780  IP=1,AXIS 
KAXIS( IP)=ieLANK 
2780  CONTINUE 

OC  2810  1=1, AXIS, 10 
KAXIS( I ) = IBAR 
2810  CONTINUE 

IF(N.LT.Il)  GOTO  2860 

IF((Y1.LE.0.0).AND.(0.0.LE.Y6))  KAXISd  Z ) = I COT 
KAXIS( IY)=I STAR 

2860  WRITE!  6,2865)  ( KA X I S ( J) , J = 1 , AX  I S ) 

2865  FORMAT! iSX, 115A1  ) 

2870  L=L+1 

IF!L.GT.!XGRI0-1 ) *10+2)  GOTO  2910 
2900  CONTINUE 
2910  RETURN 

CEBUG  SUBCHK 
END 
C 
C 

C ***FOUOIER  TRANSFORM  SUBROUT  INF#** 

C 

SUBROUTINE  CFFT2 ( A, B , NTQT ,N , NSP AN  , I SN) 
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